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Abstract 

Motion in the atmosphere or mantle convection are two among phenomena of natural 

■ convection induced by internal heat sources. They bifurcate from the conduction state 
I as a result of its loss of stability. In spite of their importance, due to the occurrence of 

■ variable coefficients in the nonlinear partial differential equations governing the evolution of 
the perturbations around the basic equilibrium, so far these phenomena were treated mostly 
numerically and experimentally. No rigorous study is known. In this paper we realize for the 

I first time such a linear study for the eigenvalue problem associated with those equations for a 

convection problem with an uniform internal heat source in a horizontal fluid layer bounded 
by two rigid walls. Our method uses Fourier series expansions for the unknown functions. 
Numerical results and graphs are given showing a destabilizing effect of the presence of the 
' heat source. 
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1 Introduction 



^ , Problems of convection induced by an internal heat source have been experimentally (e.g. 

5h I [5]) and numerically (e.g. [4]) studied by many researchers. The investigations concerned the 

effects of the heating and cooling rate, many investigations were performed for different condi- 
tions imposed on the lower and upper boundaries. In these cases it is difficult to compare the 
evaluations of the heat fluxes in the heated fluid layer with the ones from the classical convec- 
tion. The extent of the convection cells when the Rayleigh number increases are greater than 
the ones in the classical Rayleigh-Benard convection problem. 

These problems are very important for the study of the motion in the atmosphere, influencing 
all other processes that take place here. 

In this paper a horizontal layer of viscous incompressible fluid with constant viscosity and 
thermal conductivity coefficients i' and k is considered [6]. In this context, the heat and hydro- 
static transfer equations [6] are 

, - (1) 
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where 77 = const, is the heating rate, 9b, Pb and pB are the potential temperature, pressure and 
density in the basic state. In the fluid, the temperature at all point varies at the same rate as 
the boundary temperature, so the problem is characterized by a constant potential temperature 



difference between the lower and the upper boundaries A^g = — 6bi- Taking into account 
this leads to the following formula for the potential temperature distribution [6] 



Ml 



y) 



(3) 



In nondimensional variables the system of equations characterizing the problem is 
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-Vp' + AU + Gre'k, 
= 0, 
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where U = (u, y, z) is the velocity, 0' and p' are the temperature and pressure deviations from 
the basic state [6], Gr is the Grashof number, Pr is the Prandtl number and is a dimensionless 
parameter characterizing the heating (cooling) rate of the layer. 

The boundaries are considered rigid and ideal heat conducting, so the boundary conditions 
read 



U 



at z 



-— and z = —. 
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2 Eigenvalue problem 

Our analytical study has two major steps: we deduce the eigenvalue problem associated 
with ©-© and then, we solve this problem using a modified method based on Fourier series 
expansions of Chandrasekhar [1] type. In [6] the numerical investigations concerned the vertical 
distribution of the total heat fluxes and their individual components for small and moderate 
supercritical Rayleigh number in the presence of a uniform heat source. 

In order to deduce the eigenvalue problem let us consider the viscous incompressible fluid 
confined into a rectangular box bounded by two rigid walls: V : {) < x < ai, {) < y < a2-, 

— —< z < —. First we deduce the eigenvalue problem governing the linear stability equivalent 

to this problem. We assume that any unknown function in (0} is of the form from [3] 
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m = , n = , where ai = — , 02 = — , L and I are the box sizes. Here m > 1 and n > 1 

ai 02 H H 

are the number of cells in the x and the y direction. Using ® we have u = u{z)e^^'^^~^^y\ v = 

-f^^y(mx+ny) ^ ^ ^ -^^y(mx+ny) ^ qi ^ o(^z)e'^rnx+ny) ^ Applying to the equation Qi the curl 

operator, the system @ becomes 



{D'^v — a?Dv) — in{D'^w — a^w) — inGrO 
D'^u — a^Du — im{D'^w — o?w) — inGrO - 
imiPP' — a?yu — inijy^ — a^)u = 0, 
(L>2 _ a^)0Pr-^ + (1 - Nz)w = 0, 
imu + inv + DW = 0, 



= 0, 
0, 



(7) 



where = + and the boundary conditions are 

u = v = w = e = Oatz = ±^. (8) 

Prom Qs, taking into account the boundary conditions (jHl), we get u = — v. Replacing this 

n 
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. [ i ' 1 ■ I : in (ITjls, we obtain v = —^DW. In this way, the system Q and the boundary conditions 
((SI) reduce to the following two-point problem 



{D^ - o?fw - a^GrO = 0, 
(D^ _ a'^)e + Pr(l - Nz)w = 



if = Did = = at 2; = — and z = -. (10) 

3 Secular equation 

Taking into account the boundary conditions imposed for the normal component of the 
velocity and for the normal component of the derivative of the velocity, a method based on 
Fourier series expansions upon complete sets of orthogonal trigonometric functions that satisfy 
all boundary conditions it is not possible to apply. In this way, at first, we expanded the unknown 
function w upon a complete set of orthogonal Chandrasekhar type functions (see Chandrasekhar 
[1]), but here these functions introduced an extraperiodicity to the problem leading to the lose 
of one of the physical parameters, e.g. the heating (cooling) rate A^. 

In order to avoid this problem, let us modify the system Q by a translation of the variable 

z, X = z + such that the eigenvalue problem becomes 

' {D^ -a^fw = a^Gre, 
{D^ - 0^)9 = -Pr{Ni - Nx)w, ^ ^ 

N 

where A^i = 1 + — , with the boundary conditions 

= DIZJ = ^ = at X = and 1. (12) 

Eliminating the unknown function w between the two equations from ()11|). the following 
differential equation in the unknown function 9 only 

(L>2 - a'^fe + Ra{Ni - Nx)9 = (13) 

is obtained, with 9 = [D"^ - a^)9 = D{D^ - a^)9 = at x = and 1. 

Using the same method described in detail in [1], we introduce a new function ^ and rewrite 
the equation (|13|1 as a new system in the unknown functions 9 and ^ 



(L>2 _ a^ fo = (AT^ -Nx)^, 
^ = -a'^Rad, 



2^ -a (14) 



where Ra = Gr ■ Pr is the Rayleigh number. The Rayleigh number is an eigenvalue in H14|l 
and {9, ^) represent the corresponding eigenvector. We are interested in finding the smallest 



eigenvalue in (|14() . defining the neutral manifold. Taking into account the boundary conditions 

_ oo _ 

for 9, we set ^ = ^ \E'fc sin /cvrx, = X^^^ ^fcQfc) where is found from the equation 



m=l 
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(L>^ - a^yQk = (Ni - Nx) sin knx. (15) 
Equation (fT3)l has the solution 
1 



G 



>27r2 + a2)3 



(^0 + + A2x'^) cosh ax + {Bq + Sj^x + i?f x^) sinh ax+ 
+( A'l — Jyx) sm /cvrx — ,, „ „ ^ cos kttx , 



The unknown coefficients Aq, A^, A2, Bq, Bi, B2 are obtained from the system followed by 
imposing to the function Gfc to satisfy the boundary conditions imposed on 6 



ik27r2 + a2' 



(^g + + ^^) cosh a + {B^ + + B|) sinh a = ^^^^^^^ (-1) 
+ aSf = -2k-KN, 

a sinh + A2 (cosh a + 2a sinh a) + a cosh ai? J + -B^ (sinh a + 2a cosh a) = (17) 
a'^A\ + 3a5f = ^k7r{k'^TT'^ + a^), 



a2 cosh aAi + ^2 cosh a + 3a sinh a) + a2 sinh aBf + i?2 (3a cosh a + 2a2 sinh a) 
= ^YL^A:4_l)fc(A:2vr2 + a2). 



Replacing the obtained expressions for the coefficients B^, i = 0, 1,2 in the equation (|14j) 9 
and imposing the condition that the obtained equation be orthogonal to the functions sinlirz, 

I G N with respect to the inner product of ^2(0,1), (/, g) = / fgdz, we obtain an infinite 

Jo 

system of algebraic linear equations in ^'fc. 

The secular equation is obtained by imposing the condition that the determinant of this 
system to vanish, i.e. 

^'*llf = (^2^2 +"a2)3 [-^0^01 + ^l^^l + ^24l + B',42 + i^l^/J^ + ^14 + N,^f- 

(18) 
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' - if A; = / 
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where T^i = < 



2fe;[(-l)fe+' - 1] 



if 7^ 




Tr{l-k){l + k) 



I' 



X* cosh(ax) sin Inxdx, if j = 1 



/ X* sinh(ax) sin lirxdx, if j = 2, i = 0, 1, 2. 
. Jo 



From the secular equation, keeping some of the parameters fixed, we obtain the neutral mani- 
folds, e.g neutral curves or neutral surfaces. 



4 Numerical results 

In this case, a first approximation for the Rayleigh number Ra is obtained by taking k = I = 1 
and imposing to the first-order minor in the matrix associated with the infinite algebraic system 
in to vanish [1] . But for A; = Z = 1 the expression of the Rayleigh number following from the 
secular equation is given by 

(7r^ + a^)5(a + sinh(a)) 

a2[(sinh(a) -|- a){a'^ + ir'^f - 8air'^{l + cosh(a))] ' 

In this expression the parameter N does not occur so, we cannot use this approximation. We 
performed numerical evaluations for the Rayleigh number for the second and the third approx- 
imation, i.e. k = I = 2 and k = I = 3 respectively. The numerical results are presented in Table 
1. 

From the approximate numerical evaluations it can be seen that when the wavenumber 
a is increasing, the approximate Rayleigh number is also increasing. When the wavenumber 
is kept constant an increase in the heating (cooling) rate parameter leads to a decreasing of 
the Rayleigh number. Moreover, for values of the parameter N greater than a certain value 
(stabilized around 10), the third approximation is greater, yet not very different, from the one 
obtained in the second approximation. 

In fig.l the approximate neutral curve for the parameter characterizing the heating rate 
= 1 is presented. When both parameters are varying the approximate neutral surface is 
obtained (fig. 2). 
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9.711 


1715.079324 


1715.079324 
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9.711 


1711.742588 
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9.711 


1701.891001 
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10.0 


1712.257687 
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10.0 


1664.341789 


1659.087870 


1701.62704 
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12.0 


1685.422373 


1680.142966 


1723.62407 
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12.0 


1547.460446 


1546.437526 


1590.19681 


9 


12.0 


1508.147637 


1508.070640 


1551.72378 


10 


12.0 


1468.449223 


1469.223193 


1512.69203 


12 


12 


1389.837162 


1392.166660 


1434.90396 
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12 


1243.442054 


1247.400439 


1288.50149 


10 


9.0 


1482.527042 


1482.391699 


1525.59302 


11 


9.0 


1446.915467 


1447.488638 


1490.55802 


12 


9.00 


1411.401914 


1412.610226 


1455.48233 


30 


12.00 


879.104231 


884.196861 


917.507873 


50 


9.00 


643.4478727 


647.561188 


673.848081 



Table 1 Numerical evaluations of the Rayleigh number for various values of the parameters 

N and a. 




2 2,5 3 3,5 4 4,5 5 



Figure 1: The approximate neutral surface Ra = Ra{a) for N = 1. 

We also verified the numerical results with a variational method: the Rayleigh quotient 
method. The obtained numerical results are also presented in Table 1. However the approxi- 
mations are limited by the difficult evaluation of the associated matrix for a large number of 
functions in the expansion sets. Also, in this cases, the convergence it is not always assured. 
The advantage in using this method is that the expressions of the neutral manifolds are easy to 
obtain. 

When N = the problem reduces to the particular case of Rayleigh-Benard convection and 
the numerical evaluation lead us to a critical value similar to the classical value for the Rayleigh 
number, i.e. Ra = 1715.079324 for a = 3.117. 



Figure 2: The approximate neutral surface Ra = Ra{a,N). 



5 Conclusions 

In this paper wc arc concerned with the stabihty of the mechanical equilibrium of a fluid 
layer, namely a problem of convection with uniform internal heat source in a parallelipipedic 
box. The eigenvalue problem is obtained and then investigated in order to obtain the secular 
equation which gives the approximate stability limits. 

Numerical investigations were performed in order to establish the influence of the two pa- 
rameters, the wavenumber and the parameter characterizing the heating (cooling) rate on the 
approximate values of the Rayleigh number: when the heating (cooling) rate is increasing, the 
domain of stability is decreasing and an increase in the wavenumber enlarge the domain of sta- 
bility. We pointed out that the first approximation does not work in this case and remarked 
that some other classical method do not work either in this case. 
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